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ANALYSIS OF TWO-EQUATION TURBULENCE MODELS 
FOR RECIRCULATING FLOWS 


S. Thangam^ 

ICASE, NASA Langley Research Center 
Hampton, Virginia 23665 


ABSTRACT 

The two-equation K-z model is used to analyze turbulent separated flow past a backward-fac- 
ing step. It is shown that if the model constants are modified to be consistent with the accepted 
energy decay rate for isotropic turbulence, the dominant features of the flow field — namely — the 
size of the separation bubble and the streamwise component of the mean velocity, can be accurately 
predicted. In addition, except in the vicinity of the step, very good predictions for the turbulent 
shear stress, the wall pressure and the wall shear stress are obtained. The model is also shown to 
provide good predictions for the turbulence intensity in the region downstream of the reattachment 
point. Estimated long-time growth rates for the turbulent kinetic energy and dissipation rate of ho- 
mogeneous shear flow are utilized to develop an optimal set of constants for the two equation K-z 
model. The physical implications of the model performance are also discussed. 
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1. INTRODUCTION 


The accurate prediction of turbulent separated flows is crucial to the analysis of many physical 
systems. One of the most commonly used approaches is to Reynolds average the equations of mo- 
tion and apply a one-point Reynolds stress closure (Launder & Spalding, 1974; Speziale, 1991). 
Among the various methods used for one-point closure, the two-equation turbulence models are 
probably the most popular since they offer an acceptable compromise between the more accurate 
but computationally expensive full Reynolds stress closure models and the less rigorous one-equa- 
tion or algebraic models. 

In its standard form the two-equation K-e model consists of a representation for the eddy vis- 
cosity in terms of the turbulent kinetic energy and the turbulent dissipation rate which are them- 
selves represented in terms of modeled transport equations (Launder & Spalding, 1974). The K-e 
model involves the introduction of five constants which include the proportionality constant in the 
eddy viscosity, followed by the two constants in the equation for the turbulent dissipation, and two 
additional constants which represent the ratios of the eddy viscosity to the diffusivities of the tur- 
bulent kinetic energy and the turbulent dissipation. Typically these constants are calibrated based 
on the available experimental findings for some simplified flow configurations (Launder & Spal- 
ding, 1974; Rodi, 1980). 

A widely used benchmark test for studying the performance of two-equation turbulence models 
involves the physical configuration of an abrupt expansion in a channel — the backward-facing 
step (cf., figure 1). The flow separates at the comer and is characterized by the presence of a large 
recirculation region which is straddled by a shear layer. The separated flow reattaches at a down- 
stream location x r and is followed by a flow recovery region. For fully-developed turbulent flow, 
an attached boundary layer exists in the region adjacent to the upper wall. It is the presence of such 
diverse features that had prompted many researchers in the past to use this flow configuration as a 
benchmark test case for analyzing the predictive capability of turbulence models. In particular, a 
variety of two-equation turbulence models have been tested and compared with the experimental 
data of Kim, Kline & Johnston (1980) and Eaton & Johnston (1981) and a description of these and 
other related results may be found in Kline, Cantwell & Johnston (1981). 
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It is also generally believed that the standard K-z model (Launder & Spalding, 1974), with wall 
functions, underpredicts the reattachment point by a substantial amount of the order of 20-25%. 
To overcome this deficiency several alternative forms of the K-z model have been developed over 
the years. Among these, Sindir (1984) made modifications to account for streamline curvature 
based on the algebraic stress model of Gibson (1985) and obtained a modest amount of improve- 
ment. Hanjalic, Launder, & Schiestel (1981) proposed a multiple scale K-z model wherein the 
turbulent kinetic energy K and the turbulent dissipation rate e were decomposed into low and high 
wavenumber parts and such a model was used by Chen (1986) to obtain significantly improved re- 
sults for the separated flow past a backward-facing step. The multiple scale representation is also 
the crux of the models based on the renormalization group theory (Yakhot & Orszag, 1986). On 
the other hand Speziale & Ngo (1989) reported comparable improvements for the backward-fac- 
ing step problem based on an anisotropic K-z model implying that the main source of the errors 
could be due to the use of an isotropic eddy-viscosity in the standard K-z model. However, Avva, 
Kline & Ferziger (1988) have suggested that the large underprediction of the reattachment point 
attributed to the standard K-z model was mainly due to the under-resolution of the computational 
domain. 

The present work is primarily aimed at the development of modifications for the standard K-z 
model to improve its predictive capability. It is shown that if the standard K-z model is modified 
to properly represent the decay rate of turbulent kinetic energy in isotropic turbulence its predictive 
capability is considerably enhanced. Computations are performed based on a finite-volume meth- 
od and it is demonstrated that the model can accurately predict the dominant features of the flow 
field. These include the size of the separation bubble, the velocity profiles, the wall pressure and 
the wall shear stress — quantities which are of considerable use from the engineering point of view. 
In this context, the criteria for an optimal choice of model constants based on the growth rates of 
turbulent kinetic energy and dissipation is developed. The physical implications of these findings 
are also discussed. 
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2. FORMULATION OF THE PHYSICAL PROBLEM 


The problem to be considered is the fully-developed turbulent flow of an incompressible vis- 
cous fluid past a backward-facing step. A schematic of the flow field is shown in figure 1 . The 
incoming flow separates in the vicinity of the step comer and reattaches at a distance x r . The inlet 
channel has a length L i and a height hi while the channel downstream of the step has a length L c 
and a height /ig. The governing equations consist of the Reynolds averaged Navier-Stokes and the 
continuity equations which may be expressed as: 


Bu Bv 
Bx + By 


Bu - Bu Bp ( B 2 u 9 2 m"\ dx dx 

- + u^- + v^- = --L-+v — - + — - 1-— 

Bt Bx By Bx (^2 3^,2 J Bx By 


Bv Bp ( B 2 v B 2 v\ 

Tt +U te +V Bj = ~Bj + \y x 2 + tf) 




( 1 ) 

( 2 ) 

(3) 


where, u and v are the mean velocity components in the x and y directions; p is the modified mean 
pressure; x , x and x are the components of the Reynolds stress tensor x. . = and v is the 

xx xy yy ij i j 

kinematic viscosity. In the standard K—e model with isotropic eddy- viscosity the Reynolds stress 
tensor takes the form (see, Launder and Spalding, 1974; Rodi, 1980) 


9 C”2 _ 

x.. = ±-K5..-2C a — S if 
ij 3 ij n e lJ 


(4) 


where 


s «-i 
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1 Bx i + Bx- 

v J 1 J 


(5) 


is the mean rate of strain tensor, K=-x u is the turbulent kinetic energy, e is the turbulent dissi- 


pation rate, and = (u,v) is the mean velocity vector. The governing equations for K and £ 
may be modeled by the following transport equations (Launder & Spalding, 1974). 
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where. 


is the eddy viscosity. 



fP= -t 


du ,du dv^ 

t ( H ) 

xx dx x y dy dx 



( 8 ) 


(9) 


is the turbulence production, and C ^ C eh C e2 , a K and a e are dimensionless constants. At this 
point a brief discussion of the method by which these constants are obtained is in order. The quan- 
tities and a e are the ratios of the eddy viscosity to the diffusivity of the turbulent kinetic energy 
and turbulent dissipation, respectively. Both quantities are of order 1 and for the standard model 
the values of a K and a e are estimated to be about 1 and 1.3, respectively, based on the results from 
two-dimensional shear flows. To evaluate C e 2 3 simple model based on homogeneous grid gener- 
ated turbulence (where production and diffusion are absent) is used. Based on the available exper- 
imental evidence and the results from direct simulations a value between 1.8 and 2.0 is recom- 
mended (with 1.92 being the recommended value for the standard K-e model). The quantity C^, 
is typically evaluated based on two-dimensional shear layers that are in local equilibrium. For this 
case, production and dissipation balance each other and from (6), (8) & (9) it can be shown that 

= -Uv/K. The measurements of Bradshaw, Ferriss & Atwell (1975) have shown that 

uv/K~ -0.3 for such a case and hence a value of = 0.09 is typically used. Next, C el is eval- 
uated from (7) by assuming that in the near- wall regions production approximately equals dissipa- 
tion and the convection of dissipation is negligible. Under these conditions with the added assump- 
tion of the logarithmic velocity profile it can be shown that 
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^el ~ ^e2 


( 10 ) 


K 2 

a tJC~v 


Herein the von Kdrmdn constant, k = 0.41, and the above expression is used to estimate ~ 1 .44. 

The constants for the standard K-e model will now be reconsidered in an effort to improve its 
predictive capability. In this regard it should be noted that the first term in the transport equation 
for turbulent dissipation (7) is the production of dissipation while the second term represents the 
destruction of dissipation. The constant C e 2 which appears in the destruction term also plays a cru- 
cial role in the decay of turbulence. For the case of isotropic turbulence the decay of turbulent ki- 
netic energy can be shown to be (Reynolds, 1987) 


K(t) = K a 


-,-l/(C e2 -l) 


l+(C e 2 -l)^ 


K 


( 11 ) 


which is consistent with the experimentally observed decay rate of K~ t 1-2 (Comte-Bellot & 
Corrsin, 1971) when C e 2 = 1 1/6. The standard K-e model on the other hand, uses C e 2 = 1.92 which 

corresponds to a power-law decay with K~t~ 11 . In addition, more recent studies on the modeling 
of pressure-strain correlations for the full Reynolds stress closure recommend the asymptotically 

—1 9 

correct decay rate of K~t ' (which again corresponds to C E 2 = 1 1/6) for homogeneous shear 
flows (Speziale, Sarkar & Gatski, 1991). Thus, C e 2 is set to 1 1/6 in the modified K-e model while 
all other constants are left unchanged. 

For the above model, the Reynolds averaged equations (l)-(3), (6) and (7) are to be solved sub- 
ject to the following boundary conditions: 

(a) inlet profiles for u, K and e are specified five step heights upstream of the step comer. For this 
purpose, a separate Reynolds- averaged calculation for the developing turbulent channel flow 
is performed. The results from these computations at an appropriate location near the outlet of 
the channel (determined by matching with the experimental data of Eaton & Johnston, 1980) 
are used as the conditions at the inlet. 

(b) Conservative extrapolated outflow conditions are applied thirty step heights downstream of the 
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step comer; these conditions involve the following: i) the u-component of the velocity for the 
cells at the outflow boundary are obtained by extrapolation; ii) the « -component of the velocity 
is then computed by the application of a mass balance; and iii) the scalar quantities such as 
pressure, turbulent kinetic energy and turbulent dissipation are all obtained by extrapolation. 
It was found that a downstream channel length of about thirty stepheights was needed to ensure 
that the local error for all the quantities was of the same order as the interior values. 

(c) At the upper and the lower walls and along the step the law of the wall is applied in the standard 
two-layer form, wherein 


3/4 K 


. 3/2 


u + = - lny + + 5, e = C, 

k J ^ y 

and the normal derivative of K is taken to vanish at the wall. These conditions are applied at 
the first grid point y away from the wall if y + =yu z /\ > 11.6 given that u + = u/u z (u^ is 


the shear velocity and k = 0.41 is the von K&mdn constant); and if y + < 11.6, then u, K and 
£ are interpolated to their wall values based on viscous sublayer constraints. It should be noted 
that the law of the wall does not formally apply to separated turbulent boundary layers. How- 
ever, since the separation point is fixed at the comer of the backstep and since the flowfield is 
solved iteratively with the shear velocity Ux updated until convergence, major errors do not ap- 
pear to result from its use (Avva et ah, 1988; Speziale & Ngo, 1989). 

The governing equations (l)-(3), (6) and (7) with the boundary conditions are discretized based 
on a finite volume method and applied for the flow past a backward-facing step. The resulting sys- 
tem of algebraic equations are solved iteratively by a line relaxation method with the repeated ap- 
plication of the tridiagonal matrix solution algorithm (see, Lilley & Rhode, 1980). Computations 
are performed for the configuration wherein the expansion ratio, £ is 1:3 (i.e., step height to outlet 
channel height, h:h 2 ) and the Reynolds number Re = 132,000 (based on the inlet centerline mean 

velocity and outlet channel height). This configuration was selected based on the fact that it has 
been used by a number of previous researchers to calibrate their turbulence models (Kline et al., 
1981) and because of the available experimental data (Kim et ah, 1980; Eaton & Johnston, 1981). 


The issue of adequate numerical resolution for this particular configuration has been consid- 
ered in the past by several researchers (Avva et ah, 1988; Thangam & Hur, 1991). Based on their 
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findings a finite- volume scheme with a 200x100 nonuniform mesh (which is known to yield results 
that are within 0.3% of the grid independent solution) was used in the present study. The particular 
type of nonuniform computational mesh used for this study is shown in figure 2 (wherein, for clar- 
ity only the alternate grid boundaries are shown over a shortened region); and for this mesh, the 

magnitude ofy + of the cells adjacent to the solid boundaries has a maximum value of about 6. The 
computed solution was assumed to have converged to its steady state when the magnitude of the 
relative average difference between successive iterations for all the velocity components, pressure, 
the mass residual, and the location of the reattachment point was less than 10' 4 . Approximately 
2000 iterations were needed for the convergence of the standard K-e model which corresponds to 
about 18 minutes of CPU time in a partially vectorized mode on the CRAY-YMP supercomputer 
using 64-bit precision. In the following section the computational results are presented and dis- 
cussed. 
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3. DISCUSSION OF THE RESULTS 


The system of governing equations outlined in §2 are solved using the modified K-e model for 
the case of flow past a backward-facing step of 1:3 expansion ratio (step height to downstream 
channel height ratio) at a flow Reynolds number of 132,000 (based on the inlet centerline velocity 
and the outlet channel height). The results of the computations are presented in figures 3-5 and 
compared with the experimental data of Kim et al. (1980) and Eaton & Johnston (1981). 

In figure 3 (a), the contours of meanflow streamlines are shown for the modified K-e model. 
As can be seen, the computed streamlines show a reattachment at XjJh ~ 7.0, which is in excellent 
agreement with the experimental mean reattachment point (~ 7.0). This is in contrast to most of 
the earlier results based on the standard K-e model which underpredict the reattachment point by 
as much as 20% (Kline et al., 1981). In fact, even the substantially more elaborate nonlinear mod- 
els have difficulty in predicting the reattachment point accurately (Speziale, 1991). In the present 
study, calculations were performed using the same computational mesh and flow conditions for the 
standard K-e model (wherein, C e2 = 1.92) to yield a reattachment point XjJh ~ 6.0. Thus, a change 
of about 5% in the value of the Ce2 (from 1.92 to 1 1/6) causes a 15% reduction in the reattachment 
length. 

We now consider the profiles of the streamwise component of the mean velocity. In figure 
3 (b), the predicted values of the streamwise mean velocity profiles are compared with the experi- 
mental data. As can be seen, there is very good agreement between the computations and the ex- 
perimental findings over the entire flow field. 

The profiles of the dimensionless turbulence shear stress are shown next in figure 4 (a) at se- 
lected locations in the streamwise direction. As can be seen the turbulence shear stresses are also 
well predicted and in good agreement with the experimental results. In figure 4 (b) the turbulence 
intensity profiles are shown and compared with the experimental results. While the model signif- 
icantly underpredicts the normal stresses in the vicinity of the reattachment point, the overall agree- 
ment can still be considered good. 

In figures 5 (a)-(b) the variation of the experimental and computed wall pressure coefficient, 
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Cp (defined as C p = 2 (p - p r ) /puf where, p r and U r are the pressure and velocity, on the cen- 
terline at the inlet) are shown along the top and the bottom walls downstream of the step for the 
modified K-e model. As can be seen, there is very good agreement between the experimental and 
computational findings. 

Another parameter of importance is the wall shear stress which could be expressed in its di- 
mensionless form as = 2 xj p U^. Since, the experimental data for the particular configuration 

used is not available in the literature, the scaled data of Driver & Seegmiller (1985) which corre- 
sponds to an expansion ratio of 1 :9 is used. It is known that except in the region very close to the 
step wall the variation of the dimensionless wall shear stress ratio, Cf/ Cf oui (wherein, Cf out corre- 
sponds to the fully-developed value) with the normalized distance, ( x-x r )/x r is essentially indepen- 
dent of the expansion ratio (Adams & Eaton, 1988). Based on this premise, the experimental re- 
sults (with an expansion ratio of 1:9 and a reattachment point, x r / h ~ 6.26) are scaled and com- 
pared with computational results (with an expansion ratio of 1:3 and a reattachment point, x r l h ~ 
7.0) and shown in figure 5 (c). As can be seen the computational results are in good agreement 
with the experimental findings. 

We now consider some of the physical implications of the results presented so far. It should 
be noted at the outset that the two-equation K-e model used here is based on the assumption that 
the eddy viscosity is isotropic (i.e., same for all components of the Reynolds stress). Thus the re- 
sults obtained clearly show that for even for such diversified flows as the separated flow past the 
backward -facing step which is characterized by a large recirculating region (wherein the normal- 
stresses and shear-stresses are of the same order, and where both are considerably smaller than the 
inertial quantities) the isotropic turbulence model is adequate. In fact, most of the physical quan- 
tities of interest could be accurately computed based on the modified K-e model, although the 
anisotropic K-e models or the full Reynolds stress closure model should be preferred for analyzing 
the recirculating region, particularly in the vicinity of the step (Speziale, 1991). 

Since the flow past the backward-facing step is characterized by a shear layer that straddles a 
large recirculation zone, the following interpretation for the size of the separation bubble is con- 
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sidered next. In a homogeneous shear flow, equations (6) and (7) can be combined to yield a long- 
time behavior for the turbulence kinetic energy K and dissipation e of the form (Speziale & Mac 
Giolla Mhuiris, 1989) 


K~exp ( Xt *) 
e~exp (X£*) 

where t* is the time (nondimensionalized by the shear rate) and X is the growth rate given by 


X = 


C„(C e2 -C el ) 2 ■ 
(C.,-1) (C,,-l) 


1/2 


( 12 ) 


Therefore, the eddy viscosity also has the form 

\ T ~exp (Xt*) 

Hence the growth rate of the eddy viscosity is intimately tied to the three model constants, C C el 
and C e 2 for shear flows. Purely for pedagogical reasons as well as for consistency (such as equa- 
tion 10) the value of C e i is maintained between 1.40 and 1.45, although there are several two equa- 
tion models wherein C E i can have a significantly different value (Yakhot & Orszag, 1986). For a 
specified value of C e i an increase in or C e 2 would cause an increase in the growth rate X, and 

therefore in Vj. An increase in v T would make the flow field more dissipative and would be ex- 
pected to reduce the size of the separation bubble. 

In this regard, it should be noted that for the standard K-e model (with = 0.09, C e j = 1.44, 
C e 2 = 1 .92), the growth rate X ~ 0.225 and the reattachment point x r /h ~ 6.0; whereas for the mod- 
ified K-e model (with = 0.09, C e \ = 1.44, C e2 = 1 1/6), the corresponding growth rate X ~ 0.195 
with the reattachment point XjJh ~ 7.0. Calculations for other values of C e2 indicate the same 
trend.; for example, when C e 2 = 1.70 (corresponding to the growth rate X ~ 0.140), x r /h ~ 9.65. 
Thus, the separation bubble size is nearly in inverse proportion to the growth rate. In addition, 
several additional computations were performed for different values of C^, C e iand C e 2 such that 
for each set the growth rate X ~ 0.2. For all these cases, the separation bubble size XjJh ~ 7.0 and 
the mean velocity profiles were observed to be in agreement with the experimental values. How- 
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ever, the profiles of the turbulence stresses as well as intensity tend to be significantly affected 
when the model constants deviate too far (i.e., > 15%) from that of the standard K-z model. In this 
context it should be pointed out that the magnitude of a# and a £ themselves do not significantly 

(i.e., of the order of 1% in terms of the mean flow quantities) affect the above findings so long as 
they are not appreciably different from unity. In the present study, this premise was verified by 
performing computations with several values of a# and a e ranging from 0.7 to 1.4. 

These findings however, should not be taken to imply as an unequivocal endorsement for the 
standard K-z model. The two equation models are based on the assumption that the local charac- 
teristics of turbulence can be represented by a single velocity scale and that the individual compo- 
nents of the Reynolds stress are related to this scale through an eddy-viscosity relationship. Such 
an approach leads to a model that cannot account for the multiple scales present in more complex 
flows. Furthermore, there are certain flow situations where a particular aspect of the motion may 
be solely due to the anisotropy of the Reynolds stress components (for example, the turbulent sec- 
ondary motion in straight, noncircular ducts), and for such flows it is important to model accurately 
the anisotropy of the Reynolds stress (Speziale, 1991). In addition, the model is purely dissipative 
and is unable to account for the relaxational effects of the Reynolds stress. Partly to overcome such 
shortcomings models based on anisotropic generalization of the eddy viscosity have been devel- 
oped (for an overview, see, Speziale; 1991) and successfully applied for the prediction of turbulent 
separated flows (Speziale & Ngo, 1989). In general such models tend to be complex and in some 
instances have been known to be numerically dispersive to the extent of requiring special treatment 
(Speziale & Ngo, 1989). 

In this regard, it should be noted that very close to the step wall another, much weaker, sec- 
ondary recirculating layer is present (see, figure 3). The simple isotropic eddy-viscosity model out- 
lined here cannot be expected to predict accurately the details, shape, or the size of such features 
of the flow field. Consequently, the pressure coefficients, the wall shear stress and the shape of the 
streamlines would be affected in the vicinity of the step. However, what this study has shown is 
that in spite of these shortcomings, the simple two-equation model based on isotropic eddy viscos- 
ity with optimally selected model constants can accurately predict the dominant features of recir- 
culating flow past the backward-facing step. 
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4. CONCLUSIONS 


An analysis of two-equation turbulence models of the K-e type is presented to improve then- 
predictive capability for separated flows. The constants of the K-e model are optimally selected 
based on the accepted energy decay rate in isotropic turbulence and then applied to the prediction 
of turbulent separated flow past the backward-facing step. 

The model is shown to accurately predict the dominant features of the flow field. The well- 
predicted features include the size of the separation bubble and the mean velocity. In addition, ex- 
cept in the vicinity of the step, the wall pressure and the wall shear stress are also well predicted. 
Furthermore, the model is shown to provide very good predictions for the turbulence shear stress, 
while its predictive capability for the turbulence intensity is good except in the region near the re- 
attachment point. 

Since the turbulent separated flow past a backward-facing step is characterized by the presence 
of a large recirculation region that is straddled by a shear layer, an empirically derived constraint 
for the model constants is proposed based on the observed long-time growth rate of turbulent ki- 
netic energy and dissipation rate in homogeneous shear layers. 

In summary, the findings of this study indicate that with minor modifications the standard two- 
equation K-e turbulence model based on isotropic eddy viscosity can be a viable option for the pre- 
diction of turbulent separated flows past the backward-facing step. 
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